# erin baggott carter
# october 14, 2019
# replication file for tables and figures

library(stargazer)
library(ggplot2)
library(systemfit)
library(tseries)

figdir = '~/Dropbox/Diversionary Aggression Paper/II Final Submission/Figures/'
setwd('~/Dropbox/Diversionary Aggression Paper/II Final Submission/Replication Files')
dat = read.csv('FINALDATA.csv',stringsAsFactors=F);head(dat)


#####################################
#        STATIONARITY TESTS         #
#####################################

# table 1 (levels stationary)
adf.test(na.omit(dat$approve))
adf.test(na.omit(dat$Lvcon))
adf.test(na.omit(dat$LvconEnemy))

# table 2 (levels not all stationary, first differences stationary)
adf.test(na.omit(dat$demApprove)) 
adf.test(na.omit(dat$indApprove)) 
adf.test(na.omit(dat$repApprove)) #

adf.test(na.omit(dat$demApprove[dat$dem==1])) # 
adf.test(na.omit(dat$indApprove[dat$dem==1])) 
adf.test(na.omit(dat$repApprove[dat$dem==1]))  

adf.test(na.omit(dat$demApprove[dat$rep==1])) #
adf.test(na.omit(dat$indApprove[dat$rep==1])) 
adf.test(na.omit(dat$repApprove[dat$rep==1]))  

adf.test(na.omit(dat$DdemApprove)) 
adf.test(na.omit(dat$DindApprove))
adf.test(na.omit(dat$DrepApprove)) 

adf.test(na.omit(dat$DdemApprove[dat$dem==1])) 
adf.test(na.omit(dat$DindApprove[dat$dem==1])) 
adf.test(na.omit(dat$DrepApprove[dat$dem==1]))  

adf.test(na.omit(dat$DdemApprove[dat$rep==1])) 
adf.test(na.omit(dat$DindApprove[dat$rep==1])) 
adf.test(na.omit(dat$DrepApprove[dat$rep==1])) 

# table 3 (levels stationary)
adf.test(na.omit(dat$mcon))
adf.test(na.omit(dat$unrate))
adf.test(na.omit(dat$misery)) 

######################################################################
# TABLE 1: HOSTILE FOREIGN POLICY RHETORIC AND PRESIDENTIAL APPROVAL #
######################################################################

m1 = lm(approve ~ Lvcon + as.factor(administration) ,data=dat); summary(m1)
m2 = lm(approve ~ Lapprove + Lvcon + Lmcon  + Lmisery + Lsecondterm  + Lunited,data=dat); summary(m2)
m3 = lm(approve ~ Lapprove + Lvcon + Lmcon  + Lmisery + Lsecondterm  + Lunited +factor(administration),data=dat); summary(m3)
m4 = lm(approve ~ LvconEnemy + as.factor(administration),data=dat); summary(m4)
m5 = lm(approve ~ Lapprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited,data=dat); summary(m5)
m6 = lm(approve ~ Lapprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited  +factor(administration),data=dat); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Hostile Foreign Policy Rhetoric and Presidential Approval",label="Approval",
          order=c("Lvcon","LvconEnemy","Lapprove","Lmcon","Lmisery","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Approval$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Misery index$_{t-1}$', 'Second term$_{t-1}$', 'Unified government$_{t-1}$'),
          dep.var.labels=c("Approval$_{t}$"),font.size="small",
          omit="administration")


#######################################################################
# FIGURE 5: HOSTILE FOREIGN POLICY RHETORIC AND PRESIDENTIAL APPROVAL #
#######################################################################

# simulate non-fixed effects model:

newdata2 <- data.frame(LvconEnemy = rep(seq(from = min(na.omit(dat$LvconEnemy)), to = max(na.omit(dat$LvconEnemy)), length.out = 7)))

newdata2$Lapprove 	= mean(na.omit(dat$Lapprove))
newdata2$Lmcon 		= mean(na.omit(dat$Lmcon))
newdata2$Lmisery 	= mean(na.omit(dat$Lmisery))
newdata2$Lsecondterm = mean(na.omit(dat$Lsecondterm))
newdata2$Lunited 	= mean(na.omit(dat$Lunited))

newdata2 = cbind(newdata2,predict(m5,newdata2,type="response",se=TRUE)) # pred probability

newdata2 = within(newdata2, {
  PredApprove = fit # exp(fit)  # change from glm -- no exp
  LL = fit - 1.96*se.fit #exp(fit - 1.96*se.fit)
  UL = fit + 1.96*se.fit #exp(fit + 1.96*se.fit)
})			

ggplot(newdata2, aes(LvconEnemy, PredApprove)) +
  geom_ribbon(aes(ymin=LL,ymax=UL),alpha=.1) + 
  geom_line(aes(),size=2) + 
  labs(x = "Rhetorical conflict in t-1", y = "Predicted approval in t")	 + 
  theme_bw() + 
  geom_vline(xintercept=c(mean(na.omit(dat$vconEnemy))), linetype="dotted") + 
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18)) 

dev.print(file=paste(figdir,"approvalPlot4.pdf",sep=""), device=pdf); dev.off()


#########################################################
# ONLINE APPENDIX FIGURE 1: RHETORIC BY ADMINISTRATION  #
#########################################################

### RHETORICAL CONFLICT

sub = setNames(data.frame(matrix(ncol = 1+length(unique(dat$administration)), nrow = nrow(dat))), c("ym",as.character(as.vector(unique(dat$administration)))))
sub$ym = dat$ym
dat$administration=as.character(as.vector(dat$administration))

vars = unique(dat$administration)
for(v in vars){
  print(v)
  for(i in 1:nrow(sub)){
    value = dat$vcon[dat$ym==sub$ym[i] & dat$administration==v]
    if(length(value)!=0){
      sub[,v][i] = value
    }
  }
}

#show 50% confidence regions and color each variable separately
sub$ym=NULL
error.bars(sub,alpha=.90,
           main="Presidential Rhetorical Conflict",col=rainbow(ncol(sub)),
           xlab=NA,ylab="Hostile Rhetoric Per Month",
           las=2,bty="n"
           )

dev.print(file=paste(figdir,"violinVcon.pdf",sep=""), device=pdf); dev.off()


### RHETORICAL COOPERATION

sub = setNames(data.frame(matrix(ncol = 1+length(unique(dat$administration)), nrow = nrow(dat))), c("ym",as.character(as.vector(unique(dat$administration)))))
sub$ym = dat$ym
dat$administration=as.character(as.vector(dat$administration))

vars = unique(dat$administration)
for(v in vars){
  print(v)
  for(i in 1:nrow(sub)){
    value = dat$vcoop[dat$ym==sub$ym[i] & dat$administration==v]
    if(length(value)!=0){
      sub[,v][i] = value
    }
  }
}

#show 50% confidence regions and color each variable separately
sub$ym=NULL
error.bars(sub,alpha=.90,
           main="Presidential Rhetorical Cooperation",col=rainbow(ncol(sub)),
           xlab=NA,ylab="Cooperative Rhetoric Per Month",
           las=2,bty="n"
)
dev.print(file=paste(figdir,"violinVcoop.pdf",sep=""), device=pdf); dev.off()

##################################################################
# TABLE 2: HOSTILE FOREIGN POLICY RHETORIC AND PARTISAN APPROVAL #
##################################################################

# democratic presidents
m1 = lm(demApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m1)
m2 = lm(demApprove ~ LdemApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1,]); summary(m2)
m3 = lm(indApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m3)
m4 = lm(indApprove ~ LindApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1,]); summary(m4)
m5 = lm(repApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m5)
m6 = lm(repApprove ~ LrepApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$dem==1,]); summary(m6)

# republican presidents

m7 = lm(repApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m7)
m8 = lm(repApprove ~ LrepApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$rep==1,]); summary(m8)
m9 = lm(indApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m9)
m10 = lm(indApprove ~ LindApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1,]); summary(m10)
m11 = lm(demApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m11)
m12 = lm(demApprove ~ LdemApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1,]); summary(m12)


stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Hostile Foreign Policy Rhetoric and Partisan Approval",
          label="partisanApproval",font.size="tiny",
          omit="administration",
          order=c("LvconEnemy","LdemApprove","LindApprove","LrepApprove","Lmcon","Lmisery","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'D approval$_{t-1}$',
                             'I approval$_{t-1}$',
                             'R approval$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Misery index$_{t-1}$', 
                             'Second term$_{t-1}$', 
                             'Unified government$_{t-1}$'))

#############################################################
# TABLE 3: MISERY INDEX AND HOSTILE FOREIGN POLICY RHETORIC #
#############################################################

# all presidents (aggregate masks differential partisan effects)
m1 = lm(vcon ~ Lmisery + as.factor(administration), data=dat); summary(m1)
m2 = lm(vcon ~ Lvcon + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m2)
m3 = lm(vconEnemy ~ Lmisery + as.factor(administration), data=dat); summary(m3)
m4 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m4)
m5 = lm(mcon ~ Lmisery + as.factor(administration), data=dat); summary(m5)
m6 = lm(mcon ~ Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Misery Index and Hostile Foreign Policy Rhetoric",label="diversion",
          dep.var.labels=c("Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$'),
          order=c("Lmisery","Lvcon","LvconEnemy","Lunited","Lmcon","Lsecondterm"),
          covariate.labels=c('Misery index$_{t-1}$',
                             'Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material conflict$_{t-1}$',
                              'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="footnotesize",
          omit="administration")

#### dems only
m1 = lm(vcon ~ Lmisery + as.factor(administration), data=dat[dat$dem==1,]); summary(m1)
m2 = lm(vcon ~ Lvcon + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$dem==1,]); summary(m2)
m3 = lm(vconEnemy ~ Lmisery + as.factor(administration), data=dat[dat$dem==1,]); summary(m3)
m4 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$dem==1,]); summary(m4)
m5 = lm(mcon ~ Lmisery + as.factor(administration), data=dat[dat$dem==1,]); summary(m5)
m6 = lm(mcon ~ Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$dem==1,]); summary(m6)

### reps only
m7 = lm(vcon ~ Lmisery + as.factor(administration), data=dat[dat$rep==1,]); summary(m7)
m8 = lm(vcon ~ Lvcon + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$rep==1,]); summary(m8)
m9 = lm(vconEnemy ~ Lmisery + as.factor(administration), data=dat[dat$rep==1,]); summary(m9)
m10 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$rep==1,]); summary(m10)
m11 = lm(mcon ~ Lmisery + as.factor(administration), data=dat[dat$rep==1,]); summary(m11)
m12 = lm(mcon ~ Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$rep==1,]); summary(m12)

stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Misery Index and Hostile Foreign Policy Rhetoric",label="diversion",
          dep.var.labels=c("Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$',"Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$'),
          order=c("Lmisery","Lvcon","LvconEnemy","Lunited","Lmcon","Lsecondterm"),
          covariate.labels=c('Misery index$_{t-1}$',
                             'Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="tiny",
          omit="administration")


##############################################################
# FIGURE 6: MISERY INDEX AND HOSTILE FOREIGN POLICY RHETORIC #
##############################################################

# simulate based on model 4

dems = as.factor(c("Truman","Kennedy","Johnson","Carter","Clinton","Obama"))	

### democratic presidents

newdata2 <- data.frame(
  #L2unrate = rep(seq(from = 2.5, to = 10.5, by =.5), 6),
  Lmisery = rep(seq(from = 2.51, to = 10.267, by =.5), 6),
  administration = factor(rep(1:6, each = length(rep(seq(from = 2.51, to = 10.267, by =.5)))), levels = 1:6, labels = levels(dems)))  # choose here

newdata2$LvconEnemy = mean(na.omit(dat$LvconEnemy))
newdata2$Lmcon 		= mean(na.omit(dat$Lmcon))
newdata2$Lsecondterm = mean(na.omit(dat$Lsecondterm))
newdata2$Lelection 	= mean(na.omit(dat$Lelection))
newdata2$Lunited 	= mean(na.omit(dat$Lunited))

newdata2 = cbind(newdata2,predict(m4,newdata2,type="response",se=TRUE)) # pred probability

newdata2 = within(newdata2, {
  PredvconEnemy = fit # exp(fit)  # change from glm -- no exp
  LL = fit - 1.96*se.fit #exp(fit - 1.96*se.fit)
  UL = fit + 1.96*se.fit #exp(fit + 1.96*se.fit)
})			


ggplot(newdata2, aes(Lmisery, PredvconEnemy)) +
  geom_ribbon(aes(ymin=LL,ymax=UL,fill=administration),alpha=.1) + 
  geom_line(aes(colour = administration),size=2) + 
  labs(x = "Misery index in t-1", y = "Hostile rhetoric toward rivals in t")	 + 
  theme_bw() + ylim(0,3.5) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18)) +
  theme(legend.title=element_blank(),
        legend.text=element_text(size=18))

dev.print(file=paste(figdir,"diversionPlotDems2.pdf",sep=""), device=pdf); dev.off()


#############################################################################
# TABLE 8: SEM MODEL OF DIVERSION, UNEMPLOYMENT, AND APPROVAL FOR DEMOCRATS #
#############################################################################

# diversion model: diversion is caused by unemployment
eq1 = vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration)
# approval model: approval is caused by unemployment and diversion
eq2 = approve ~ Lapprove + LvconEnemy + Lmisery + Lmcon + Lsecondterm  + Lunited +factor(administration)

sub = dat[c('dem','vconEnemy', 'LvconEnemy','Lmcon','Lmisery','Lsecondterm','Lelection','Lunited','administration','approve','Lapprove')]
eqSystem = list(a = eq1, b = eq2)
fitols = systemfit(eqSystem, data=na.omit(sub))
options(digits=2)
options(scipen=100)
summary(fitols)

#####################################################################################################
# ONLINE APPENDIX TABLE 2: HOSTILE FOREIGN POLICY RHETORIC AND PRESIDENTIAL APPROVAL (UNEMPLOYMENT) #
#####################################################################################################

m1 = lm(approve ~ Lvcon + as.factor(administration),data=dat); summary(m1)
m2 = lm(approve ~ Lapprove + Lvcon + Lmcon  + Lunrate + Lsecondterm  + Lunited,data=dat); summary(m2)
m3 = lm(approve ~ Lapprove + Lvcon + Lmcon  + Lunrate + Lsecondterm  + Lunited + factor(administration),data=dat); summary(m3)
m4 = lm(approve ~ LvconEnemy +factor(administration),data=dat); summary(m4)
m5 = lm(approve ~ Lapprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm  + Lunited,data=dat); summary(m5)
m6 = lm(approve ~ Lapprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm  + Lunited  +factor(administration),data=dat); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Hostile Foreign Policy Rhetoric and Presidential Approval (Unemployment)",
          label="ApprovalUnemployment",
          order=c("Lvcon","LvconEnemy","Lapprove","Lmcon","Lunrate","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Approval$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Unemployment$_{t-1}$', 'Second term$_{t-1}$', 'Unified government$_{t-1}$'),
          dep.var.labels=c("Approval$_{t}$"),font.size="small",
          omit="administration")


#####################################################################################################
# ONLINE APPENDIX TABLE 3: HOSTILE FOREIGN POLICY RHETORIC AND PARTISAN APPROVAL (UNEMPLOYMENT)     #
#####################################################################################################

# democratic presidents
m1 = lm(demApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m1)
m2 = lm(demApprove ~ LdemApprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1,]); summary(m2)
m3 = lm(indApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m3)
m4 = lm(indApprove ~ LindApprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1,]); summary(m4)
m5 = lm(repApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m5)
m6 = lm(repApprove ~ LrepApprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$dem==1,]); summary(m6)

# republican presidents
m7 = lm(repApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m7)
m8 = lm(repApprove ~ LrepApprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$rep==1,]); summary(m8)
m9 = lm(indApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m9)
m10 = lm(indApprove ~ LindApprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1,]); summary(m10)
m11 = lm(demApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m11)
m12 = lm(demApprove ~ LdemApprove + LvconEnemy + Lmcon + Lunrate + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1,]); summary(m12)


stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Hostile Foreign Policy Rhetoric and Partisan Approval (Unemployment)",
          label="partisanApprovalUnemployment",font.size="tiny",
          omit="administration",
          order=c("LvconEnemy","LdemApprove","LindApprove","LrepApprove","Lmcon","Lunrate","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'D approval$_{t-1}$',
                             'I approval$_{t-1}$',
                             'R approval$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Unemployment$_{t-1}$', 
                             'Second term$_{t-1}$', 
                             'Unified government$_{t-1}$'))

#####################################################################################################
# ONLINE APPENDIX TABLE 4: UNEMPLOYMENT AND HOSTILE FOREIGN POLICY RHETORIC                         #
#####################################################################################################

# aggregate
m1 = lm(vcon ~ Lunrate + as.factor(administration), data=dat); summary(m1)
m2 = lm(vcon ~ Lvcon + Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m2)
m3 = lm(vconEnemy ~ Lunrate + as.factor(administration), data=dat); summary(m3)
m4 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m4)
m5 = lm(mcon ~ Lunrate + as.factor(administration), data=dat); summary(m5)
m6 = lm(mcon ~ Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Unemployment and Hostile Foreign Policy Rhetoric",label="diversionUnemployment",
          dep.var.labels=c("Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$'),
          order=c("Lunrate","Lvcon","LvconEnemy","Lunited","Lmcon","Lsecondterm"),
          covariate.labels=c('Unemployment$_{t-1}$',
                             'Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="small",
          omit="administration")

# dem only
m1 = lm(vcon ~ Lunrate + as.factor(administration), data=dat[dat$dem==1,]); summary(m1)
m2 = lm(vcon ~ Lvcon + Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$dem==1,]); summary(m2)
m3 = lm(vconEnemy ~ Lunrate + as.factor(administration), data=dat[dat$dem==1,]); summary(m3)
m4 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$dem==1,]); summary(m4)
m5 = lm(mcon ~ Lunrate + as.factor(administration), data=dat[dat$dem==1,]); summary(m5)
m6 = lm(mcon ~ Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$dem==1,]); summary(m6)

# rep only
m7 = lm(vcon ~ Lunrate + as.factor(administration), data=dat[dat$rep==1,]); summary(m7)
m8 = lm(vcon ~ Lvcon + Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$rep==1,]); summary(m8)
m9 = lm(vconEnemy ~ Lunrate + as.factor(administration), data=dat[dat$rep==1,]); summary(m9)
m10 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$rep==1,]); summary(m10)
m11 = lm(mcon ~ Lunrate + as.factor(administration), data=dat[dat$rep==1,]); summary(m11)
m12 = lm(mcon ~ Lmcon + Lunrate +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$rep==1,]); summary(m12)


stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Unemployment and Hostile Foreign Policy Rhetoric",label="diversionUnemployment",
          dep.var.labels=c("Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$',"Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$'),
          order=c("Lunrate","Lvcon","LvconEnemy","Lunited","Lmcon","Lsecondterm"),
          covariate.labels=c('Unemployment$_{t-1}$',
                             'Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="tiny",
          omit="administration")


##################################################################################################
# ONLINE APPENDIX TABLE 5: HOSTILE FOREIGN POLICY RHETORIC AND PRESIDENTIAL APPROVAL (EX NIXON)  #
##################################################################################################

m1 = lm(approve ~ Lvcon + as.factor(administration) ,data=dat[dat$administration!="Nixon",]); summary(m1)
m2 = lm(approve ~ Lapprove + Lvcon + Lmcon  + Lmisery + Lsecondterm  + Lunited,data=dat[dat$administration!="Nixon",]); summary(m2)
m3 = lm(approve ~ Lapprove + Lvcon + Lmcon  + Lmisery + Lsecondterm  + Lunited +factor(administration),data=dat[dat$administration!="Nixon",]); summary(m3)
m4 = lm(approve ~ LvconEnemy + as.factor(administration),data=dat[dat$administration!="Nixon",]); summary(m4)
m5 = lm(approve ~ Lapprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited,data=dat[dat$administration!="Nixon",]); summary(m5)
m6 = lm(approve ~ Lapprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited  +factor(administration),data=dat[dat$administration!="Nixon",]); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Hostile Foreign Policy Rhetoric and Presidential Approval (Ex Nixon)",label="Approval Ex Nixon",
          order=c("Lvcon","LvconEnemy","Lapprove","Lmcon","Lmisery","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Approval$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Misery index$_{t-1}$', 'Second term$_{t-1}$', 'Unified government$_{t-1}$'),
          dep.var.labels=c("Approval$_{t}$"),font.size="small",
          omit="administration")


##################################################################################################
# ONLINE APPENDIX TABLE 6: HOSTILE FOREIGN POLICY RHETORIC AND PARTISAN APPROVAL (EX NIXON)      #
##################################################################################################

# democratic presidents
m1 = lm(demApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1 & dat$administration!="Nixon",]); summary(m1)
m2 = lm(demApprove ~ LdemApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1 & dat$administration!="Nixon",]); summary(m2)
m3 = lm(indApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1 & dat$administration!="Nixon",]); summary(m3)
m4 = lm(indApprove ~ LindApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1 & dat$administration!="Nixon",]); summary(m4)
m5 = lm(repApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$dem==1 & dat$administration!="Nixon",]); summary(m5)
m6 = lm(repApprove ~ LrepApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$dem==1 & dat$administration!="Nixon",]); summary(m6)

# republican presidents
m7 = lm(repApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1 & dat$administration!="Nixon",]); summary(m7)
m8 = lm(repApprove ~ LrepApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$rep==1 & dat$administration!="Nixon",]); summary(m8)
m9 = lm(indApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1 & dat$administration!="Nixon",]); summary(m9)
m10 = lm(indApprove ~ LindApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1 & dat$administration!="Nixon",]); summary(m10)
m11 = lm(demApprove ~ LvconEnemy + as.factor(administration),data=dat[dat$rep==1 & dat$administration!="Nixon",]); summary(m11)
m12 = lm(demApprove ~ LdemApprove + LvconEnemy + Lmcon + Lmisery + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1 & dat$administration!="Nixon",]); summary(m12)

stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Hostile Foreign Policy Rhetoric and Partisan Approval (Ex Nixon)",
          label="partisanApproval exNixon",font.size="tiny",
          omit="administration",
          order=c("LvconEnemy","LdemApprove","LindApprove","LrepApprove","Lmcon","Lmisery","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'D approval$_{t-1}$',
                             'I approval$_{t-1}$',
                             'R approval$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Misery index$_{t-1}$', 
                             'Second term$_{t-1}$', 
                             'Unified government$_{t-1}$'))


##################################################################################################
# ONLINE APPENDIX TABLE 7: MISERY INDEX AND HOSTILE FOREIGN POLICY RHETORIC (EX NIXON)           #
##################################################################################################

# aggregate
m1 = lm(vcon ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon",]); summary(m1)
m2 = lm(vcon ~ Lvcon + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon",]); summary(m2)
m3 = lm(vconEnemy ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon",]); summary(m3)
m4 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon",]); summary(m4)
m5 = lm(mcon ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon",]); summary(m5)
m6 = lm(mcon ~ Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon",]); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Misery Index and Hostile Foreign Policy Rhetoric (Ex Nixon)",label="diversion ex nixon",
          dep.var.labels=c("Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$'),
          order=c("Lmisery","Lvcon","LvconEnemy","Lunited","Lmcon","Lsecondterm"),
          covariate.labels=c('Misery index$_{t-1}$',
                             'Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="small",
          omit="administration")

# dem only
m1 = lm(vcon ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon" & dat$dem==1,]); summary(m1)
m2 = lm(vcon ~ Lvcon + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon" & dat$dem==1,]); summary(m2)
m3 = lm(vconEnemy ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon" & dat$dem==1,]); summary(m3)
m4 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon" & dat$dem==1,]); summary(m4)
m5 = lm(mcon ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon" & dat$dem==1,]); summary(m5)
m6 = lm(mcon ~ Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon" & dat$dem==1,]); summary(m6)

# rep only
m7 = lm(vcon ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon" & dat$rep==1,]); summary(m7)
m8 = lm(vcon ~ Lvcon + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon" & dat$rep==1,]); summary(m8)
m9 = lm(vconEnemy ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon" & dat$rep==1,]); summary(m9)
m10 = lm(vconEnemy ~ LvconEnemy + Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon" & dat$rep==1,]); summary(m10)
m11 = lm(mcon ~ Lmisery + as.factor(administration), data=dat[dat$administration!="Nixon" & dat$rep==1,]); summary(m11)
m12 = lm(mcon ~ Lmcon + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat[dat$administration!="Nixon" & dat$rep==1,]); summary(m12)

stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Misery Index and Hostile Foreign Policy Rhetoric (Ex Nixon)",label="diversion ex nixon",
          dep.var.labels=c("Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$',"Rhetorical conflict$_t$","Rhetorical conflict $\\rightarrow$ rivals$_t$",'Material conflict$_t$'),
          order=c("Lmisery","Lvcon","LvconEnemy","Lunited","Lmcon","Lsecondterm"),
          covariate.labels=c('Misery index$_{t-1}$',
                             'Rhetorical conflict$_{t-1}$',
                             'Rhetorical conflict$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material conflict$_{t-1}$',
                             'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="tiny",
          omit="administration")

##################################################################################################
# ONLINE APPENDIX TABLE 8: COOPERATIVE FOREIGN POLICY RHETORIC AND PRESIDENTIAL APPROVAL         #
##################################################################################################

m1 = lm(approve ~ Lvcoop + as.factor(administration) ,data=dat); summary(m1)
m2 = lm(approve ~ Lapprove + Lvcoop + Lmcoop  + Lmisery + Lsecondterm  + Lunited,data=dat); summary(m2)
m3 = lm(approve ~ Lapprove + Lvcoop + Lmcoop  + Lmisery + Lsecondterm  + Lunited +factor(administration),data=dat); summary(m3)
m4 = lm(approve ~ LvcoopEnemy + as.factor(administration),data=dat); summary(m4)
m5 = lm(approve ~ Lapprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm  + Lunited,data=dat); summary(m5)
m6 = lm(approve ~ Lapprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm  + Lunited  +factor(administration),data=dat); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Cooperative Foreign Policy Rhetoric and Presidential Approval",label="Approval Cooperative",
          order=c("Lvcoop","LvcoopEnemy","Lapprove","Lmcoop","Lmisery","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical cooperation$_{t-1}$',
                             'Rhetorical cooperation$\\rightarrow$rivals$_{t-1}$',
                             'Approval$_{t-1}$',
                             'Material cooperation$_{t-1}$',
                             'Misery index$_{t-1}$', 'Second term$_{t-1}$', 'Unified government$_{t-1}$'),
          dep.var.labels=c("Approval$_{t}$"),font.size="small",
          omit="administration")

##################################################################################################
# ONLINE APPENDIX TABLE 9: COOPERATIVE FOREIGN POLICY RHETORIC AND PARTISAN APPROVAL             #
##################################################################################################

# democratic presidents
m1 = lm(demApprove ~ LvcoopEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m1)
m2 = lm(demApprove ~ LdemApprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1,]); summary(m2)
m3 = lm(indApprove ~ LvcoopEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m3)
m4 = lm(indApprove ~ LindApprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm  + Lunited + as.factor(administration),data=dat[dat$dem==1,]); summary(m4)
m5 = lm(repApprove ~ LvcoopEnemy + as.factor(administration),data=dat[dat$dem==1,]); summary(m5)
m6 = lm(repApprove ~ LrepApprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$dem==1,]); summary(m6)

# republican presidents
m7 = lm(repApprove ~ LvcoopEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m7)
m8 = lm(repApprove ~ LrepApprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm  + Lunited  + as.factor(administration),data=dat[dat$rep==1,]); summary(m8)
m9 = lm(indApprove ~ LvcoopEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m9)
m10 = lm(indApprove ~ LindApprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1,]); summary(m10)
m11 = lm(demApprove ~ LvcoopEnemy + as.factor(administration),data=dat[dat$rep==1,]); summary(m11)
m12 = lm(demApprove ~ LdemApprove + LvcoopEnemy + Lmcoop + Lmisery + Lsecondterm + Lunited + as.factor(administration),data=dat[dat$rep==1,]); summary(m12)


stargazer(m1,m2,m3,m4,m5,m6, m7,m8,m9,m10,m11,m12,title="Cooperative Foreign Policy Rhetoric and Partisan Approval",
          label="partisanApproval cooperation",font.size="tiny",
          omit="administration",
          order=c("LvcoopEnemy","LdemApprove","LindApprove","LrepApprove","Lmcoop","Lmisery","Lsecondterm","Lunited"),
          covariate.labels=c('Rhetorical cooperation$\\rightarrow$rivals$_{t-1}$',
                             'D approval$_{t-1}$',
                             'I approval$_{t-1}$',
                             'R approval$_{t-1}$',
                             'Material cooperation$_{t-1}$',
                             'Misery index$_{t-1}$', 
                             'Second term$_{t-1}$', 
                             'Unified government$_{t-1}$'))

##################################################################################################
# ONLINE APPENDIX TABLE 10: MISERY INDEX AND COOPERATIVE FOREIGN POLICY RHETORIC                 #
##################################################################################################

m1 = lm(vcoop ~ Lmisery + as.factor(administration), data=dat); summary(m1)
m2 = lm(vcoop ~ Lvcoop + Lmcoop + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m2)
m3 = lm(vcoopEnemy ~ Lmisery + as.factor(administration), data=dat); summary(m3)
m4 = lm(vcoopEnemy ~ LvcoopEnemy + Lmcoop + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m4)
m5 = lm(mcoop ~ Lmisery + as.factor(administration), data=dat); summary(m5)
m6 = lm(mcoop ~ Lmcoop + Lmisery +Lsecondterm +Lelection  +Lunited + factor(administration), data=dat); summary(m6)

stargazer(m1,m2,m3,m4,m5,m6,title="Misery Index and Cooperative Foreign Policy",label="diversion cooperative",
          dep.var.labels=c("Rhetorical cooperation$_t$","Rhetorical cooperation $\\rightarrow$ rivals$_t$",'Material cooperation$_t$'),
          order=c("Lmisery","Lvcoop","LvcoopEnemy","Lunited","Lmcoop","Lsecondterm"),
          covariate.labels=c('Misery index$_{t-1}$',
                             'Rhetorical cooperation$_{t-1}$',
                             'Rhetorical cooperation$\\rightarrow$rivals$_{t-1}$',
                             'Unified government$_{t-1}$',
                             'Material cooperation$_{t-1}$',
                             'Second term$_{t-1}$',
                             'Election year$_{t-1}$'),
          font.size="small",
          omit="administration")


############################################################################
# FIGURE 1: NET US INTERACTIONS WITH MAIN PARTNERS, 1946-2010              #
############################################################################


par(mfrow=c(2,1))
tab = arrange(c,desc(net_mcon),country)[1:5,] # top 5 net material conflict targets (mcon - mcoop)
tab = tab[c("country","net_mcon")] # top 5 targets of material conflict 
df = data.frame(sort(as.vector(tab$net_mcon),decreasing=F)) # reverse order
rownames(df) = rev(tab$country) # reverse order
x = t(as.matrix(df))
x
barplot(x, col="coral4",ylim=c(0,200),ylab="Net material conflict events",main="Net Material Conflict Recipients, 1946-2010") 
# top 10 net verbal cooperation targets (vcoop - vcon)
tab = arrange(c,desc(net_vcoop),country)[1:5,]
tab = tab[c("country","net_vcoop")] # top 5 targets of net verbal cooperation 
df = data.frame(sort(as.vector(tab$net_vcoop),decreasing=F)) # reverse order
rownames(df) = rev(tab$country) # reverse order
x = t(as.matrix(df))
x
barplot(x, col="DarkSeaGreen",ylim=c(0,1000),ylab="Net rhetorical cooperation events",main="Net Rhetorical Cooperation Recipients, 1946-2010") 

dev.print(file=paste(figdir,"net2",sep=""), device=pdf); dev.off()


############################################################################
# FIGURE 2: US DIPLOMACY TOWARD FOUR COUNTRIES                             #
############################################################################

# load day level data
df = read.csv("AmericanDiplomacyDatasetUSinitiated.csv",stringsAsFactors=F);head(df)

########################### china

dfCHN <- as.data.frame(as.Date(as.Date("1946-01-01"):as.Date("2010-12-31"), origin="1970-01-01"))
colnames(dfCHN) <- "date"
dfCHN[c("scale","vcoop","mcoop","vcon","mcon")] <- 0

sub = df[df$target=="CHN",];nrow(sub)
for(i in 1:nrow(dfCHN)){
  if((i %% 1000==0)==T){print(paste('On',i,'of',nrow(dfCHN),'rows,',round((i/nrow(dfCHN))*100,2),'percent'))}
  dfCHN$scale[i] <- mean(sub$scale[dfCHN$date[i]==sub$date & sub$target=="CHN"])
  dfCHN$vcoop[i] <- sum(sub$vcoop[dfCHN$date[i]==sub$date & sub$target=="CHN"])
  dfCHN$mcoop[i] <- sum(sub$mcoop[dfCHN$date[i]==sub$date & sub$target=="CHN"])
  dfCHN$vcon[i] <- sum(sub$vcon[dfCHN$date[i]==sub$date & sub$target=="CHN"])
  dfCHN$mcon[i] <- sum(sub$mcon[dfCHN$date[i]==sub$date & sub$target=="CHN"])
}

annCHN = as.data.frame(c(1946:2010))
colnames(annCHN) = "year"
for(i in 1:nrow(annCHN)){
  annCHN$vcoop[i] = sum(na.omit(dfCHN$vcoop[substr(dfCHN$date,1,4)==annCHN$year[i]]))
  annCHN$mcoop[i] = sum(na.omit(dfCHN$mcoop[substr(dfCHN$date,1,4)==annCHN$year[i]]))
  annCHN$vcon[i] = sum(na.omit(dfCHN$vcon[substr(dfCHN$date,1,4)==annCHN$year[i]]))
  annCHN$mcon[i] = sum(na.omit(dfCHN$mcon[substr(dfCHN$date,1,4)==annCHN$year[i]]))
}

########################### afghanistan

dfAFG <- as.data.frame(as.Date(as.Date("1946-01-01"):as.Date("2010-12-31"), origin="1970-01-01"))
colnames(dfAFG) <- "date"
dfAFG[c("scale","vcoop","mcoop","vcon","mcon")] <- 0

sub = df[df$target=="AFG",];nrow(sub)
for(i in 1:nrow(dfAFG)){
  if((i %% 1000==0)==T){print(paste('On',i,'of',nrow(dfAFG),'rows,',round((i/nrow(dfAFG))*100,2),'percent'))}
  dfAFG$scale[i] <- mean(sub$scale[dfAFG$date[i]==sub$date & sub$target=="AFG"])
  dfAFG$vcoop[i] <- sum(sub$vcoop[dfAFG$date[i]==sub$date & sub$target=="AFG"])
  dfAFG$mcoop[i] <- sum(sub$mcoop[dfAFG$date[i]==sub$date & sub$target=="AFG"])
  dfAFG$vcon[i] <- sum(sub$vcon[dfAFG$date[i]==sub$date & sub$target=="AFG"])
  dfAFG$mcon[i] <- sum(sub$mcon[dfAFG$date[i]==sub$date & sub$target=="AFG"])
}
# make annual data
annAFG = as.data.frame(c(1946:2010))
colnames(annAFG) = "year"
for(i in 1:nrow(annAFG)){
  annAFG$vcoop[i] = sum(na.omit(dfAFG$vcoop[substr(dfAFG$date,1,4)==annAFG$year[i]]))
  annAFG$mcoop[i] = sum(na.omit(dfAFG$mcoop[substr(dfAFG$date,1,4)==annAFG$year[i]]))
  annAFG$vcon[i] = sum(na.omit(dfAFG$vcon[substr(dfAFG$date,1,4)==annAFG$year[i]]))
  annAFG$mcon[i] = sum(na.omit(dfAFG$mcon[substr(dfAFG$date,1,4)==annAFG$year[i]]))
}

########################### iraq

dfIRQ <- as.data.frame(as.Date(as.Date("1946-01-01"):as.Date("2010-12-31"), origin="1970-01-01"))
colnames(dfIRQ) <- "date"
dfIRQ[c("scale","vcoop","mcoop","vcon","mcon")] <- 0

sub = df[df$target=="IRQ",];nrow(sub)
for(i in 1:nrow(dfIRQ)){
  if((i %% 1000==0)==T){print(paste('On',i,'of',nrow(dfIRQ),'rows,',round((i/nrow(dfIRQ))*100,2),'percent'))}
  dfIRQ$scale[i] <- mean(sub$scale[dfIRQ$date[i]==sub$date & sub$target=="IRQ"])
  dfIRQ$vcoop[i] <- sum(sub$vcoop[dfIRQ$date[i]==sub$date & sub$target=="IRQ"])
  dfIRQ$mcoop[i] <- sum(sub$mcoop[dfIRQ$date[i]==sub$date & sub$target=="IRQ"])
  dfIRQ$vcon[i] <- sum(sub$vcon[dfIRQ$date[i]==sub$date & sub$target=="IRQ"])
  dfIRQ$mcon[i] <- sum(sub$mcon[dfIRQ$date[i]==sub$date & sub$target=="IRQ"])
}
# make annual data
annIRQ = as.data.frame(c(1946:2010))
colnames(annIRQ) = "year"
for(i in 1:nrow(annIRQ)){
  annIRQ$vcoop[i] = sum(na.omit(dfIRQ$vcoop[substr(dfIRQ$date,1,4)==annIRQ$year[i]]))
  annIRQ$mcoop[i] = sum(na.omit(dfIRQ$mcoop[substr(dfIRQ$date,1,4)==annIRQ$year[i]]))
  annIRQ$vcon[i] = sum(na.omit(dfIRQ$vcon[substr(dfIRQ$date,1,4)==annIRQ$year[i]]))
  annIRQ$mcon[i] = sum(na.omit(dfIRQ$mcon[substr(dfIRQ$date,1,4)==annIRQ$year[i]]))
}

########################### iran

dfIRN <- as.data.frame(as.Date(as.Date("1946-01-01"):as.Date("2010-12-31"), origin="1970-01-01"))
colnames(dfIRN) <- "date"
dfIRN[c("scale","vcoop","mcoop","vcon","mcon")] <- 0

sub = df[df$target=="IRN",];nrow(sub)
for(i in 1:nrow(dfIRN)){
  if((i %% 1000==0)==T){print(paste('On',i,'of',nrow(dfIRN),'rows,',round((i/nrow(dfIRN))*100,2),'percent'))}
  dfIRN$scale[i] <- mean(sub$scale[dfIRN$date[i]==sub$date & sub$target=="IRN"])
  dfIRN$vcoop[i] <- sum(sub$vcoop[dfIRN$date[i]==sub$date & sub$target=="IRN"])
  dfIRN$mcoop[i] <- sum(sub$mcoop[dfIRN$date[i]==sub$date & sub$target=="IRN"])
  dfIRN$vcon[i] <- sum(sub$vcon[dfIRN$date[i]==sub$date & sub$target=="IRN"])
  dfIRN$mcon[i] <- sum(sub$mcon[dfIRN$date[i]==sub$date & sub$target=="IRN"])
}
# make annual data
annIRN = as.data.frame(c(1946:2010))
colnames(annIRN) = "year"
for(i in 1:nrow(annIRN)){
  annIRN$vcoop[i] = sum(na.omit(dfIRN$vcoop[substr(dfIRN$date,1,4)==annIRN$year[i]]))
  annIRN$mcoop[i] = sum(na.omit(dfIRN$mcoop[substr(dfIRN$date,1,4)==annIRN$year[i]]))
  annIRN$vcon[i] = sum(na.omit(dfIRN$vcon[substr(dfIRN$date,1,4)==annIRN$year[i]]))
  annIRN$mcon[i] = sum(na.omit(dfIRN$mcon[substr(dfIRN$date,1,4)==annIRN$year[i]]))
}

########################### plot

# plot type of US initiated events towards 6 countries, 5 year moving average

# plot five year moving average of each event type (two sided)
#ma <- function(x,n=5){filter(x,rep(1/n,n), sides=2)} 
ma <- function(x,n=3){filter(x,rep(1/n,n), sides=2)} 


par(mfrow=c(2,2))
# RUS
#plot(annRUS$year,ma(annRUS$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,50),ylab="No. events",xlab=NA,main="Russian Federation")
#points(annRUS$year,ma(annRUS$mcoop),type="l",col="DarkSlateGrey")
#points(annRUS$year,ma(annRUS$vcon),type="l",col="coral")
#points(annRUS$year,ma(annRUS$mcon),type="l",col="coral4")
#legend("topright",pch=16,col=c('DarkSeaGreen','DarkSlateGrey','coral','coral4'),c('Verbal cooperation','Material cooperation','Verbal conflict','Material conflict'),bty="n") 
# VMN
#plot(annVMN$year,ma(annVMN$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,50),ylab="No. events",xlab=NA,main="Vietnam")
#points(annVMN$year,ma(annVMN$mcoop),type="l",col="DarkSlateGrey")
#points(annVMN$year,ma(annVMN$vcon),type="l",col="coral")
#points(annVMN$year,ma(annVMN$mcon),type="l",col="coral4")
#legend("topright",pch=16,col=c('DarkSeaGreen','DarkSlateGrey','coral','coral4'),c('Verbal cooperation','Material cooperation','Verbal conflict','Material conflict'),bty="n") 
# AFG
plot(annAFG$year,ma(annAFG$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,50),ylab="# events/yr, 3-yr 2-sided moving avg",xlab=NA,main="Afghanistan")
points(annAFG$year,ma(annAFG$mcoop),type="l",col="DarkSlateGrey")
points(annAFG$year,ma(annAFG$vcon),type="l",col="coral")
points(annAFG$year,ma(annAFG$mcon),type="l",col="coral4")
legend("topleft",pch=16,col=c('DarkSeaGreen','DarkSlateGrey','coral','coral4'),c('Rhetorical coop.','Material coop.','Rhetorical conflict','Material conflict'),bty="n",cex=.7) 
abline(v="1979",lty="dashed") # ussr invasion begins
text(x=1979,y=50,labels="USSR",cex=.7)
text(x=1979,y=47,labels="War",cex=.7)
abline(v="1989",lty="dashed") # ussr invasion ends
text(x=1989,y=50,labels="USSR",cex=.7)
text(x=1989,y=47,labels="Withdraw",cex=.7)
abline(v="2001",lty="dashed") # invasion of afghanistan
text(x=2001,y=50,labels="USA",cex=.7)
text(x=2001,y=47,labels="War",cex=.7)
# IRQ
plot(annIRQ$year,ma(annIRQ$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,50),ylab=NA,xlab=NA,main="Iraq")
points(annIRQ$year,ma(annIRQ$mcoop),type="l",col="DarkSlateGrey")
points(annIRQ$year,ma(annIRQ$vcon),type="l",col="coral")
points(annIRQ$year,ma(annIRQ$mcon),type="l",col="coral4")
abline(v="1991",lty="dashed") # first gulf war
text(x=1991,y=50,labels="Gulf",cex=.7)
text(x=1991,y=47,labels="War",cex=.7)
abline(v="2003",lty="dashed") # war in iraq
text(x=2003-3,y=50,labels="Iraq",cex=.7)
text(x=2003-3,y=47,labels="War",cex=.7)
abline(v="2007",lty="dashed") # surge
text(x=2007+1,y=50,labels="Surge",cex=.7)
# IRN
plot(annIRN$year,ma(annIRN$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,50),ylab=NA,xlab=NA,main="Iran")
points(annIRN$year,ma(annIRN$mcoop),type="l",col="DarkSlateGrey")
points(annIRN$year,ma(annIRN$vcon),type="l",col="coral")
points(annIRN$year,ma(annIRN$mcon),type="l",col="coral4")
abline(v="1979",lty="dashed") # hostage crisis begins
text(x=1979,y=50,labels="Hostage",cex=.7)
text(x=1979,y=47,labels="Crisis",cex=.7)
# CHN
plot(annCHN$year,ma(annCHN$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,50),ylab=NA,xlab=NA,main="China")
points(annCHN$year,ma(annCHN$mcoop),type="l",col="DarkSlateGrey")
points(annCHN$year,ma(annCHN$vcon),type="l",col="coral")
points(annCHN$year,ma(annCHN$mcon),type="l",col="coral4")
abline(v="1951",lty="dashed") # korean war fighting
text(x=1951-2,y=40,labels="Korean War",cex=.7,srt=90)
abline(v="1966",lty="dashed") # cultural revolution begins
text(x=1966-2,y=40,labels="Cultural Rev.",cex=.7,srt=90)
abline(v="1972",lty="dashed") # nixon visit
text(x=1972-2,y=40,labels="Nixon Visit",cex=.7,srt=90)
abline(v="1979",lty="dashed") # normalization
text(x=1979-2,y=40,labels="Normalization",cex=.7,srt=90)
abline(v="1989",lty="dashed") # tiananmen
text(x=1989-2,y=40,labels="Tiananmen",cex=.7,srt=90)
abline(v="2001",lty="dashed") # ep3
text(x=2001-2,y=40,labels="EP-3 Crisis",cex=.7,srt=90)

dev.print(file=paste(figdir,"movingAvg4Countries3YrMovingAvg.pdf",sep=""), device=pdf); dev.off()



############################################################################
# FIGURE 3: US DIPLOMACY TOWARD WORLD                                      #
############################################################################

# make annual data
ann = as.data.frame(c(1946:2010))
colnames(ann) = "year"
for(i in 1:nrow(ann)){
  ann$vcoop[i] = sum(na.omit(dat$vcoop[dat$year==ann$year[i]]))
  ann$mcoop[i] = sum(na.omit(dat$mcoop[dat$year==ann$year[i]]))
  ann$vcon[i] = sum(na.omit(dat$vcon[dat$year==ann$year[i]]))
  ann$mcon[i] = sum(na.omit(dat$mcon[dat$year==ann$year[i]]))
  #
  ann$vconAlly[i] = sum(na.omit(dat$vconAlly[dat$year==ann$year[i]]))
  ann$vconCore[i] = sum(na.omit(dat$vconCore[dat$year==ann$year[i]]))
  ann$vconGreat[i] = sum(na.omit(dat$vconGreat[dat$year==ann$year[i]]))
  ann$vconEnemy[i] = sum(na.omit(dat$vconEnemy[dat$year==ann$year[i]]))
  #
  ann$vcoopAlly[i] = sum(na.omit(dat$vcoopAlly[dat$year==ann$year[i]]))
  ann$vcoopCore[i] = sum(na.omit(dat$vcoopCore[dat$year==ann$year[i]]))
  ann$vcoopGreat[i] = sum(na.omit(dat$vcoopGreat[dat$year==ann$year[i]]))
  ann$vcoopEnemy[i] = sum(na.omit(dat$vcoopEnemy[dat$year==ann$year[i]]))
  #
  ann$vconPol[i] = sum(na.omit(dat$vconPol[dat$year==ann$year[i]]))
  ann$vconEcon[i] = sum(na.omit(dat$vconEcon[dat$year==ann$year[i]]))
  ann$mconPol[i] = sum(na.omit(dat$mconPol[dat$year==ann$year[i]]))
  ann$mconEcon[i] = sum(na.omit(dat$mconEcon[dat$year==ann$year[i]]))
  #
  ann$vcoopPol[i] = sum(na.omit(dat$vcoopPol[dat$year==ann$year[i]]))
  ann$vcoopEcon[i] = sum(na.omit(dat$vcoopEcon[dat$year==ann$year[i]]))
  ann$mcoopPol[i] = sum(na.omit(dat$mcoopPol[dat$year==ann$year[i]]))
  ann$mcoopEcon[i] = sum(na.omit(dat$mcoopEcon[dat$year==ann$year[i]]))	
}


# plot five year moving average of each event type (two sided)
ma <- function(x,n=5){filter(x,rep(1/n,n), sides=2)} 
#Types of US Initiated Events, 5 Year Moving Average
plot(ann$year,ma(ann$vcoop),type="l",col="DarkSeaGreen",ylim=c(0,500),ylab="Annual number of events, 5 year 2-sided moving average",xlab=NA)
points(ann$year,ma(ann$mcoop),type="l",col="DarkSlateGrey")
points(ann$year,ma(ann$vcon),type="l",col="coral")
points(ann$year,ma(ann$mcon),type="l",col="coral4")
legend("topright",pch=16,col=c('DarkSeaGreen','DarkSlateGrey','coral','coral4'),c('Rhetorical cooperation','Material cooperation','Rhetorical conflict','Material conflict'),bty="n") 
dev.print(file=paste(figdir,"movingAvg.pdf",sep=""), device=pdf); dev.off()

############################################################################
# FIGURE 4: US DIPLOMACY TOWARD WORLD (BY TOPIC)                           #
############################################################################

par(mfrow=c(2,1))
# political
plot(ann$year,ma(ann$vcoopPol),type="l",col="DarkSeaGreen",ylim=c(0,500),ylab="# events/yr, 5-yr 2-sided moving avg", xlab=NA, main="Political-Military Affairs")
points(ann$year,ma(ann$mcoopPol),type="l",col="DarkSlateGrey")
points(ann$year,ma(ann$vconPol),type="l",col="coral")
points(ann$year,ma(ann$mconPol),type="l",col="coral4")
legend("topright",pch=16,col=c('DarkSeaGreen','DarkSlateGrey','coral','coral4'),c('Rhetorical cooperation','Material cooperation','Rhetorical conflict','Material conflict'),bty="n") 
# economic
plot(ann$year,ma(ann$vcoopEcon),type="l",col="DarkSeaGreen",ylim=c(0,500),ylab="# events/yr, 5-yr 2-sided moving avg",xlab=NA,main="Economic Affairs")
points(ann$year,ma(ann$mcoopEcon),type="l",col="DarkSlateGrey")
points(ann$year,ma(ann$vconEcon),type="l",col="coral")
points(ann$year,ma(ann$mconEcon),type="l",col="coral4")
legend("topright",pch=16,col=c('DarkSeaGreen','DarkSlateGrey','coral','coral4'),c('Rhetorical cooperation','Material cooperation','Rhetorical conflict','Material conflict'),bty="n") 

dev.print(file=paste(figdir,"movingAvgByTopic.pdf",sep=""), device=pdf); dev.off()


###################################################################################################
# FIGURE 7: NUMBER OF GLOBAL MATERIAL CONFLICTS ACCORDING TO DIFFERENT CORPORA (COUNTS AND MEANS) # 
###################################################################################################


notopdup = read.csv('globalDiplomaticEventsNoTopicsInclDuplicates.csv');head(notopdup)
notopdup$source=substr(notopdup$source,1,3)
notopdup = notopdup[notopdup$source=="USA",]

notop = read.csv('globalDiplomaticEventsNoTopics.csv');head(notop)
notop$source=substr(notop$source,1,3)
notop = notop[notop$source=="USA",]

top = read.csv('globalDiplomaticEvents.csv');head(top)
top$source=substr(top$source,1,3)
top = top[top$source=="USA",]

notopdup$year = as.integer(substr(notopdup$date,1,4))
notop$year = as.integer(substr(notop$date,1,4))
top$year = as.integer(substr(top$date,1,4))


################################################################

# apply cameo scores (the ones that appear in this dataset)

sort(unique(top$code))
top$scale[top$code=="14O"] <- -6.5  # fix one botched code here before numericizing
top$code <- as.numeric(as.vector(top$code))
top$scale[top$code==010] <- 0
top$scale[top$code==011] <- -.1
top$scale[top$code==012] <- -.4
top$scale[top$code==013] <- .4
top$scale[top$code==014] <- 0
top$scale[top$code==015] <- 0
top$scale[top$code==016] <- 3.4
top$scale[top$code==017] <- 0
top$scale[top$code==018] <- 3.4
top$scale[top$code==019] <- 3.4#
top$scale[top$code==020] <- 3
top$scale[top$code==021] <- 3.4
top$scale[top$code==0211] <- 3.4
top$scale[top$code==0213] <- 3.4#
top$scale[top$code==0214] <- 3.4#
top$scale[top$code==022] <- 3.4
top$scale[top$code==023] <- 3.4
top$scale[top$code==0231] <- 3.4
top$scale[top$code==0232] <- 3.4
top$scale[top$code==0233] <- 3.4
top$scale[top$code==0234] <- 3.4
top$scale[top$code==024] <- -.3
top$scale[top$code==0241] <- -.3
top$scale[top$code==0242] <- -.3
top$scale[top$code==0243] <- -.3
top$scale[top$code==0244] <- -.3
top$scale[top$code==025] <- -.3
top$scale[top$code==0253] <- -.3#
top$scale[top$code==0254] <- -.3#
top$scale[top$code==0256] <- -.3#
top$scale[top$code==026] <- 4
top$scale[top$code==027] <- 4
top$scale[top$code==028] <- 4
top$scale[top$code==030] <- 4
top$scale[top$code==031] <- 5.2
top$scale[top$code==0311] <- 5.2
top$scale[top$code==0312] <- 5.2
top$scale[top$code==032] <- 4.5
top$scale[top$code==033] <- 5.2
top$scale[top$code==0331] <- 5.2
top$scale[top$code==0332] <- 5.2
top$scale[top$code==0333] <- 5.2
top$scale[top$code==0334] <- 6
top$scale[top$code==034] <- 7
top$scale[top$code==0341] <- 7
top$scale[top$code==0344] <- 7
top$scale[top$code==035] <- 7
top$scale[top$code==0351] <- 7
top$scale[top$code==0353] <- 7
top$scale[top$code==0355] <- 7
top$scale[top$code==0356] <- 7
top$scale[top$code==036] <- 4
top$scale[top$code==037] <- 5
top$scale[top$code==038] <- 7
top$scale[top$code==039] <- 5
top$scale[top$code==040] <- 1
top$scale[top$code==041] <- 1
top$scale[top$code==042] <- 1.9
top$scale[top$code==043] <- 2.8
top$scale[top$code==044] <- 2.5
top$scale[top$code==045] <- 5
top$scale[top$code==046] <- 7
top$scale[top$code==050] <- 3.5
top$scale[top$code==051] <- 3.4
top$scale[top$code==052] <- 3.5
top$scale[top$code==053] <- 3.8
top$scale[top$code==054] <- 6
top$scale[top$code==055] <- 7
top$scale[top$code==056] <- 7
top$scale[top$code==057] <- 8
top$scale[top$code==060] <- 6
top$scale[top$code==061] <- 6.4
top$scale[top$code==062] <- 7.4
top$scale[top$code==063] <- 7.4
top$scale[top$code==064] <- 7
top$scale[top$code==070] <- 7
top$scale[top$code==071] <- 7.4
top$scale[top$code==072] <- 8.3
top$scale[top$code==073] <- 7.4
top$scale[top$code==074] <- 8.5
top$scale[top$code==075] <- 7
top$scale[top$code==080] <- 5
top$scale[top$code==081] <- 5
top$scale[top$code==0811] <- 5
top$scale[top$code==0812] <- 5
top$scale[top$code==0813] <- 5
top$scale[top$code==0814] <- 5
top$scale[top$code==082] <- 5
top$scale[top$code==083] <- 5
top$scale[top$code==0831] <- 5
top$scale[top$code==0833] <- 5
top$scale[top$code==0834] <- 5
top$scale[top$code==084] <- 7
top$scale[top$code==0841] <- 7
top$scale[top$code==0842] <- 7
top$scale[top$code==085] <- 7
top$scale[top$code==086] <- 9
top$scale[top$code==0861] <- 9
top$scale[top$code==0862] <- 9
top$scale[top$code==0863] <- 9
top$scale[top$code==087] <- 9
top$scale[top$code==0871] <- 9
top$scale[top$code==0872] <- 9
top$scale[top$code==0873] <- 9
top$scale[top$code==0874] <- 10
top$scale[top$code==090] <- -2
top$scale[top$code==091] <- -2
top$scale[top$code==092] <- -2
top$scale[top$code==093] <- -2
top$scale[top$code==094] <- -2
top$scale[top$code==100] <- -5
top$scale[top$code==101] <- -5
top$scale[top$code==1013] <- -5#
top$scale[top$code==1014] <- -5#
top$scale[top$code==102] <- -5
top$scale[top$code==1031] <- -5#
top$scale[top$code==1033] <- -5#
top$scale[top$code==1041] <- -5
top$scale[top$code==1042] <- -5
top$scale[top$code==1043] <- -5
top$scale[top$code==1044] <- -5
top$scale[top$code==105] <- -5
top$scale[top$code==1053] <- -5#
top$scale[top$code==1056] <- -5#
top$scale[top$code==106] <- -5
top$scale[top$code==107] <- -5
top$scale[top$code==109] <- -5#
top$scale[top$code==110] <- -2
top$scale[top$code==111] <- -2
top$scale[top$code==112] <- -2
top$scale[top$code==1121] <- -2
top$scale[top$code==1122] <- -2
top$scale[top$code==1123] <- -2
top$scale[top$code==1124] <- -2
top$scale[top$code==1125] <- -2
top$scale[top$code==113] <- -2
top$scale[top$code==114] <- -2
top$scale[top$code==115] <- -2
top$scale[top$code==116] <- -2#
top$scale[top$code==120] <- -4
top$scale[top$code==121] <- -4
top$scale[top$code==1211] <- -4
top$scale[top$code==1212] <- -4
top$scale[top$code==1213] <- -4
top$scale[top$code==122] <- -4
top$scale[top$code==1222] <- -4#
top$scale[top$code==123] <- -4
top$scale[top$code==1231] <- -4
top$scale[top$code==1233] <- -4
top$scale[top$code==124] <- -5
top$scale[top$code==1241] <- -5#
top$scale[top$code==1243] <- -5#
top$scale[top$code==1244] <- -5#
top$scale[top$code==1246] <- -5#
top$scale[top$code==125] <- -5
top$scale[top$code==127] <- -5
top$scale[top$code==128] <- -5
top$scale[top$code==129] <- -5#
top$scale[top$code==130] <- -4.4
top$scale[top$code==131] <- -5.8
top$scale[top$code==1311] <- -5.8
top$scale[top$code==1312] <- -5.8
top$scale[top$code==1313] <- -5.8
top$scale[top$code==132] <- -5.8
top$scale[top$code==1322] <- -5.8
top$scale[top$code==133] <- -5.8
top$scale[top$code==134] <- -5.8
top$scale[top$code==136] <- -7
top$scale[top$code==137] <- -7
top$scale[top$code==138] <- -7
top$scale[top$code==1383] <- -7
top$scale[top$code==1384] <- -7
top$scale[top$code==139] <- -7
top$scale[top$code==140] <- -6.5
top$scale[top$code==141] <- -6.5
top$scale[top$code==1411] <- -6.5
top$scale[top$code==1412] <- -6.5
top$scale[top$code==1413] <- -6.5
top$scale[top$code==1414] <- -6.5
top$scale[top$code==142] <- -6.5
top$scale[top$code==143] <- -6.5
top$scale[top$code==144] <- -7.5
top$scale[top$code==145] <- -7.5
top$scale[top$code==150] <- -7.2
top$scale[top$code==151] <- -7.2
top$scale[top$code==152] <- -7.2
top$scale[top$code==153] <- -7.2
top$scale[top$code==154] <- -7.2
top$scale[top$code==160] <- -4
top$scale[top$code==161] <- -4
top$scale[top$code==162] <- -5.6
top$scale[top$code==1621] <- -5.6
top$scale[top$code==1622] <- -5.6
top$scale[top$code==163] <- -6.5
top$scale[top$code==164] <- -7
top$scale[top$code==166] <- -8
top$scale[top$code==1662] <- -8#
top$scale[top$code==170] <- -7
top$scale[top$code==171] <- -9.2
top$scale[top$code==1711] <- -9.2
top$scale[top$code==1712] <- -9.2
top$scale[top$code==172] <- -5
top$scale[top$code==1721] <- -5
top$scale[top$code==1722] <- -5
top$scale[top$code==1723] <- -5
top$scale[top$code==1724] <- -5
top$scale[top$code==173] <- -5
top$scale[top$code==174] <- -5
top$scale[top$code==175] <- -9
top$scale[top$code==180] <- -9
top$scale[top$code==181] <- -9
top$scale[top$code==182] <- -9.5
top$scale[top$code==1821] <- -9
top$scale[top$code==1822] <- -9
top$scale[top$code==1823] <- -10
top$scale[top$code==183] <- -10
top$scale[top$code==1831] <- -10
top$scale[top$code==1832] <- -10
top$scale[top$code==185] <- -8
top$scale[top$code==186] <- -10
top$scale[top$code==190] <- -10
top$scale[top$code==191] <- -9.5
top$scale[top$code==192] <- -9.5
top$scale[top$code==193] <- -10
top$scale[top$code==194] <- -10
top$scale[top$code==195] <- -10
top$scale[top$code==196] <- -9.5
top$scale[top$code==201] <- -9.5
top$scale[top$code==202] <- -10
top$scale[top$code==203] <- -10

sort(unique(top$code)) # no NAs


################################################################

# assign verbal and material cooperation labels

top$vcoop <- 0
top$mcoop <- 0
top$vcon <- 0
top$mcon <- 0
top$vcoop[top$code >=1 & top$code <=5] <- 1
top$mcoop[top$code >=6 & top$code <=9] <- 1
top$vcoop[top$code >=10 & top$code <=59] <- 1
top$mcoop[top$code >=60 & top$code <=99] <- 1
top$vcon[top$code >=100 & top$code <= 142] <- 1
top$mcon[top$code >=143 & top$code <=204] <- 1
top$vcoop[top$code >=211 & top$code <=356] <- 1
top$mcoop[top$code >=811 & top$code <= 874] <- 1
top$vcon[top$code >= 1011 & top$code <=1424] <- 1
top$mcon[top$code >=1431 & top$code <=2042] <- 1

################################################################

# apply cameo scores (the ones that appear in this dataset)

sort(unique(notop$code))
notop$scale[notop$code=="14O"] <- -6.5  # fix one botched code here before numericizing
notop$code <- as.numeric(as.vector(notop$code))
notop$scale[notop$code==010] <- 0
notop$scale[notop$code==011] <- -.1
notop$scale[notop$code==012] <- -.4
notop$scale[notop$code==013] <- .4
notop$scale[notop$code==014] <- 0
notop$scale[notop$code==015] <- 0
notop$scale[notop$code==016] <- 3.4
notop$scale[notop$code==017] <- 0
notop$scale[notop$code==018] <- 3.4
notop$scale[notop$code==019] <- 3.4#
notop$scale[notop$code==020] <- 3
notop$scale[notop$code==021] <- 3.4
notop$scale[notop$code==0211] <- 3.4
notop$scale[notop$code==0213] <- 3.4#
notop$scale[notop$code==0214] <- 3.4#
notop$scale[notop$code==022] <- 3.4
notop$scale[notop$code==023] <- 3.4
notop$scale[notop$code==0231] <- 3.4
notop$scale[notop$code==0232] <- 3.4
notop$scale[notop$code==0233] <- 3.4
notop$scale[notop$code==0234] <- 3.4
notop$scale[notop$code==024] <- -.3
notop$scale[notop$code==0241] <- -.3
notop$scale[notop$code==0242] <- -.3
notop$scale[notop$code==0243] <- -.3
notop$scale[notop$code==0244] <- -.3
notop$scale[notop$code==025] <- -.3
notop$scale[notop$code==0253] <- -.3#
notop$scale[notop$code==0254] <- -.3#
notop$scale[notop$code==0256] <- -.3#
notop$scale[notop$code==026] <- 4
notop$scale[notop$code==027] <- 4
notop$scale[notop$code==028] <- 4
notop$scale[notop$code==030] <- 4
notop$scale[notop$code==031] <- 5.2
notop$scale[notop$code==0311] <- 5.2
notop$scale[notop$code==0312] <- 5.2
notop$scale[notop$code==032] <- 4.5
notop$scale[notop$code==033] <- 5.2
notop$scale[notop$code==0331] <- 5.2
notop$scale[notop$code==0332] <- 5.2
notop$scale[notop$code==0333] <- 5.2
notop$scale[notop$code==0334] <- 6
notop$scale[notop$code==034] <- 7
notop$scale[notop$code==0341] <- 7
notop$scale[notop$code==0344] <- 7
notop$scale[notop$code==035] <- 7
notop$scale[notop$code==0351] <- 7
notop$scale[notop$code==0353] <- 7
notop$scale[notop$code==0355] <- 7
notop$scale[notop$code==0356] <- 7
notop$scale[notop$code==036] <- 4
notop$scale[notop$code==037] <- 5
notop$scale[notop$code==038] <- 7
notop$scale[notop$code==039] <- 5
notop$scale[notop$code==040] <- 1
notop$scale[notop$code==041] <- 1
notop$scale[notop$code==042] <- 1.9
notop$scale[notop$code==043] <- 2.8
notop$scale[notop$code==044] <- 2.5
notop$scale[notop$code==045] <- 5
notop$scale[notop$code==046] <- 7
notop$scale[notop$code==050] <- 3.5
notop$scale[notop$code==051] <- 3.4
notop$scale[notop$code==052] <- 3.5
notop$scale[notop$code==053] <- 3.8
notop$scale[notop$code==054] <- 6
notop$scale[notop$code==055] <- 7
notop$scale[notop$code==056] <- 7
notop$scale[notop$code==057] <- 8
notop$scale[notop$code==060] <- 6
notop$scale[notop$code==061] <- 6.4
notop$scale[notop$code==062] <- 7.4
notop$scale[notop$code==063] <- 7.4
notop$scale[notop$code==064] <- 7
notop$scale[notop$code==070] <- 7
notop$scale[notop$code==071] <- 7.4
notop$scale[notop$code==072] <- 8.3
notop$scale[notop$code==073] <- 7.4
notop$scale[notop$code==074] <- 8.5
notop$scale[notop$code==075] <- 7
notop$scale[notop$code==080] <- 5
notop$scale[notop$code==081] <- 5
notop$scale[notop$code==0811] <- 5
notop$scale[notop$code==0812] <- 5
notop$scale[notop$code==0813] <- 5
notop$scale[notop$code==0814] <- 5
notop$scale[notop$code==082] <- 5
notop$scale[notop$code==083] <- 5
notop$scale[notop$code==0831] <- 5
notop$scale[notop$code==0833] <- 5
notop$scale[notop$code==0834] <- 5
notop$scale[notop$code==084] <- 7
notop$scale[notop$code==0841] <- 7
notop$scale[notop$code==0842] <- 7
notop$scale[notop$code==085] <- 7
notop$scale[notop$code==086] <- 9
notop$scale[notop$code==0861] <- 9
notop$scale[notop$code==0862] <- 9
notop$scale[notop$code==0863] <- 9
notop$scale[notop$code==087] <- 9
notop$scale[notop$code==0871] <- 9
notop$scale[notop$code==0872] <- 9
notop$scale[notop$code==0873] <- 9
notop$scale[notop$code==0874] <- 10
notop$scale[notop$code==090] <- -2
notop$scale[notop$code==091] <- -2
notop$scale[notop$code==092] <- -2
notop$scale[notop$code==093] <- -2
notop$scale[notop$code==094] <- -2
notop$scale[notop$code==100] <- -5
notop$scale[notop$code==101] <- -5
notop$scale[notop$code==1013] <- -5#
notop$scale[notop$code==1014] <- -5#
notop$scale[notop$code==102] <- -5
notop$scale[notop$code==1031] <- -5#
notop$scale[notop$code==1033] <- -5#
notop$scale[notop$code==1041] <- -5
notop$scale[notop$code==1042] <- -5
notop$scale[notop$code==1043] <- -5
notop$scale[notop$code==1044] <- -5
notop$scale[notop$code==105] <- -5
notop$scale[notop$code==1053] <- -5#
notop$scale[notop$code==1056] <- -5#
notop$scale[notop$code==106] <- -5
notop$scale[notop$code==107] <- -5
notop$scale[notop$code==109] <- -5#
notop$scale[notop$code==110] <- -2
notop$scale[notop$code==111] <- -2
notop$scale[notop$code==112] <- -2
notop$scale[notop$code==1121] <- -2
notop$scale[notop$code==1122] <- -2
notop$scale[notop$code==1123] <- -2
notop$scale[notop$code==1124] <- -2
notop$scale[notop$code==1125] <- -2
notop$scale[notop$code==113] <- -2
notop$scale[notop$code==114] <- -2
notop$scale[notop$code==115] <- -2
notop$scale[notop$code==116] <- -2#
notop$scale[notop$code==120] <- -4
notop$scale[notop$code==121] <- -4
notop$scale[notop$code==1211] <- -4
notop$scale[notop$code==1212] <- -4
notop$scale[notop$code==1213] <- -4
notop$scale[notop$code==122] <- -4
notop$scale[notop$code==1222] <- -4#
notop$scale[notop$code==123] <- -4
notop$scale[notop$code==1231] <- -4
notop$scale[notop$code==1233] <- -4
notop$scale[notop$code==124] <- -5
notop$scale[notop$code==1241] <- -5#
notop$scale[notop$code==1243] <- -5#
notop$scale[notop$code==1244] <- -5#
notop$scale[notop$code==1246] <- -5#
notop$scale[notop$code==125] <- -5
notop$scale[notop$code==127] <- -5
notop$scale[notop$code==128] <- -5
notop$scale[notop$code==129] <- -5#
notop$scale[notop$code==130] <- -4.4
notop$scale[notop$code==131] <- -5.8
notop$scale[notop$code==1311] <- -5.8
notop$scale[notop$code==1312] <- -5.8
notop$scale[notop$code==1313] <- -5.8
notop$scale[notop$code==132] <- -5.8
notop$scale[notop$code==1322] <- -5.8
notop$scale[notop$code==133] <- -5.8
notop$scale[notop$code==134] <- -5.8
notop$scale[notop$code==136] <- -7
notop$scale[notop$code==137] <- -7
notop$scale[notop$code==138] <- -7
notop$scale[notop$code==1383] <- -7
notop$scale[notop$code==1384] <- -7
notop$scale[notop$code==139] <- -7
notop$scale[notop$code==140] <- -6.5
notop$scale[notop$code==141] <- -6.5
notop$scale[notop$code==1411] <- -6.5
notop$scale[notop$code==1412] <- -6.5
notop$scale[notop$code==1413] <- -6.5
notop$scale[notop$code==1414] <- -6.5
notop$scale[notop$code==142] <- -6.5
notop$scale[notop$code==143] <- -6.5
notop$scale[notop$code==144] <- -7.5
notop$scale[notop$code==145] <- -7.5
notop$scale[notop$code==150] <- -7.2
notop$scale[notop$code==151] <- -7.2
notop$scale[notop$code==152] <- -7.2
notop$scale[notop$code==153] <- -7.2
notop$scale[notop$code==154] <- -7.2
notop$scale[notop$code==160] <- -4
notop$scale[notop$code==161] <- -4
notop$scale[notop$code==162] <- -5.6
notop$scale[notop$code==1621] <- -5.6
notop$scale[notop$code==1622] <- -5.6
notop$scale[notop$code==163] <- -6.5
notop$scale[notop$code==164] <- -7
notop$scale[notop$code==166] <- -8
notop$scale[notop$code==1662] <- -8#
notop$scale[notop$code==170] <- -7
notop$scale[notop$code==171] <- -9.2
notop$scale[notop$code==1711] <- -9.2
notop$scale[notop$code==1712] <- -9.2
notop$scale[notop$code==172] <- -5
notop$scale[notop$code==1721] <- -5
notop$scale[notop$code==1722] <- -5
notop$scale[notop$code==1723] <- -5
notop$scale[notop$code==1724] <- -5
notop$scale[notop$code==173] <- -5
notop$scale[notop$code==174] <- -5
notop$scale[notop$code==175] <- -9
notop$scale[notop$code==180] <- -9
notop$scale[notop$code==181] <- -9
notop$scale[notop$code==182] <- -9.5
notop$scale[notop$code==1821] <- -9
notop$scale[notop$code==1822] <- -9
notop$scale[notop$code==1823] <- -10
notop$scale[notop$code==183] <- -10
notop$scale[notop$code==1831] <- -10
notop$scale[notop$code==1832] <- -10
notop$scale[notop$code==185] <- -8
notop$scale[notop$code==186] <- -10
notop$scale[notop$code==190] <- -10
notop$scale[notop$code==191] <- -9.5
notop$scale[notop$code==192] <- -9.5
notop$scale[notop$code==193] <- -10
notop$scale[notop$code==194] <- -10
notop$scale[notop$code==195] <- -10
notop$scale[notop$code==196] <- -9.5
notop$scale[notop$code==201] <- -9.5
notop$scale[notop$code==202] <- -10
notop$scale[notop$code==203] <- -10

sort(unique(notop$code)) # no NAs


################################################################

# assign verbal and material cooperation labels

notop$vcoop <- 0
notop$mcoop <- 0
notop$vcon <- 0
notop$mcon <- 0
notop$vcoop[notop$code >=1 & notop$code <=5] <- 1
notop$mcoop[notop$code >=6 & notop$code <=9] <- 1
notop$vcoop[notop$code >=10 & notop$code <=59] <- 1
notop$mcoop[notop$code >=60 & notop$code <=99] <- 1
notop$vcon[notop$code >=100 & notop$code <= 142] <- 1
notop$mcon[notop$code >=143 & notop$code <=204] <- 1
notop$vcoop[notop$code >=211 & notop$code <=356] <- 1
notop$mcoop[notop$code >=811 & notop$code <= 874] <- 1
notop$vcon[notop$code >= 1011 & notop$code <=1424] <- 1
notop$mcon[notop$code >=1431 & notop$code <=2042] <- 1

################################################################

################################################################

# apply cameo scores (the ones that appear in this dataset)

sort(unique(notopdup$code))
notopdup$scale[notopdup$code=="14O"] <- -6.5  # fix one botched code here before numericizing
notopdup$code <- as.numeric(as.vector(notopdup$code))
notopdup$scale[notopdup$code==010] <- 0
notopdup$scale[notopdup$code==011] <- -.1
notopdup$scale[notopdup$code==012] <- -.4
notopdup$scale[notopdup$code==013] <- .4
notopdup$scale[notopdup$code==014] <- 0
notopdup$scale[notopdup$code==015] <- 0
notopdup$scale[notopdup$code==016] <- 3.4
notopdup$scale[notopdup$code==017] <- 0
notopdup$scale[notopdup$code==018] <- 3.4
notopdup$scale[notopdup$code==019] <- 3.4#
notopdup$scale[notopdup$code==020] <- 3
notopdup$scale[notopdup$code==021] <- 3.4
notopdup$scale[notopdup$code==0211] <- 3.4
notopdup$scale[notopdup$code==0213] <- 3.4#
notopdup$scale[notopdup$code==0214] <- 3.4#
notopdup$scale[notopdup$code==022] <- 3.4
notopdup$scale[notopdup$code==023] <- 3.4
notopdup$scale[notopdup$code==0231] <- 3.4
notopdup$scale[notopdup$code==0232] <- 3.4
notopdup$scale[notopdup$code==0233] <- 3.4
notopdup$scale[notopdup$code==0234] <- 3.4
notopdup$scale[notopdup$code==024] <- -.3
notopdup$scale[notopdup$code==0241] <- -.3
notopdup$scale[notopdup$code==0242] <- -.3
notopdup$scale[notopdup$code==0243] <- -.3
notopdup$scale[notopdup$code==0244] <- -.3
notopdup$scale[notopdup$code==025] <- -.3
notopdup$scale[notopdup$code==0253] <- -.3#
notopdup$scale[notopdup$code==0254] <- -.3#
notopdup$scale[notopdup$code==0256] <- -.3#
notopdup$scale[notopdup$code==026] <- 4
notopdup$scale[notopdup$code==027] <- 4
notopdup$scale[notopdup$code==028] <- 4
notopdup$scale[notopdup$code==030] <- 4
notopdup$scale[notopdup$code==031] <- 5.2
notopdup$scale[notopdup$code==0311] <- 5.2
notopdup$scale[notopdup$code==0312] <- 5.2
notopdup$scale[notopdup$code==032] <- 4.5
notopdup$scale[notopdup$code==033] <- 5.2
notopdup$scale[notopdup$code==0331] <- 5.2
notopdup$scale[notopdup$code==0332] <- 5.2
notopdup$scale[notopdup$code==0333] <- 5.2
notopdup$scale[notopdup$code==0334] <- 6
notopdup$scale[notopdup$code==034] <- 7
notopdup$scale[notopdup$code==0341] <- 7
notopdup$scale[notopdup$code==0344] <- 7
notopdup$scale[notopdup$code==035] <- 7
notopdup$scale[notopdup$code==0351] <- 7
notopdup$scale[notopdup$code==0353] <- 7
notopdup$scale[notopdup$code==0355] <- 7
notopdup$scale[notopdup$code==0356] <- 7
notopdup$scale[notopdup$code==036] <- 4
notopdup$scale[notopdup$code==037] <- 5
notopdup$scale[notopdup$code==038] <- 7
notopdup$scale[notopdup$code==039] <- 5
notopdup$scale[notopdup$code==040] <- 1
notopdup$scale[notopdup$code==041] <- 1
notopdup$scale[notopdup$code==042] <- 1.9
notopdup$scale[notopdup$code==043] <- 2.8
notopdup$scale[notopdup$code==044] <- 2.5
notopdup$scale[notopdup$code==045] <- 5
notopdup$scale[notopdup$code==046] <- 7
notopdup$scale[notopdup$code==050] <- 3.5
notopdup$scale[notopdup$code==051] <- 3.4
notopdup$scale[notopdup$code==052] <- 3.5
notopdup$scale[notopdup$code==053] <- 3.8
notopdup$scale[notopdup$code==054] <- 6
notopdup$scale[notopdup$code==055] <- 7
notopdup$scale[notopdup$code==056] <- 7
notopdup$scale[notopdup$code==057] <- 8
notopdup$scale[notopdup$code==060] <- 6
notopdup$scale[notopdup$code==061] <- 6.4
notopdup$scale[notopdup$code==062] <- 7.4
notopdup$scale[notopdup$code==063] <- 7.4
notopdup$scale[notopdup$code==064] <- 7
notopdup$scale[notopdup$code==070] <- 7
notopdup$scale[notopdup$code==071] <- 7.4
notopdup$scale[notopdup$code==072] <- 8.3
notopdup$scale[notopdup$code==073] <- 7.4
notopdup$scale[notopdup$code==074] <- 8.5
notopdup$scale[notopdup$code==075] <- 7
notopdup$scale[notopdup$code==080] <- 5
notopdup$scale[notopdup$code==081] <- 5
notopdup$scale[notopdup$code==0811] <- 5
notopdup$scale[notopdup$code==0812] <- 5
notopdup$scale[notopdup$code==0813] <- 5
notopdup$scale[notopdup$code==0814] <- 5
notopdup$scale[notopdup$code==082] <- 5
notopdup$scale[notopdup$code==083] <- 5
notopdup$scale[notopdup$code==0831] <- 5
notopdup$scale[notopdup$code==0833] <- 5
notopdup$scale[notopdup$code==0834] <- 5
notopdup$scale[notopdup$code==084] <- 7
notopdup$scale[notopdup$code==0841] <- 7
notopdup$scale[notopdup$code==0842] <- 7
notopdup$scale[notopdup$code==085] <- 7
notopdup$scale[notopdup$code==086] <- 9
notopdup$scale[notopdup$code==0861] <- 9
notopdup$scale[notopdup$code==0862] <- 9
notopdup$scale[notopdup$code==0863] <- 9
notopdup$scale[notopdup$code==087] <- 9
notopdup$scale[notopdup$code==0871] <- 9
notopdup$scale[notopdup$code==0872] <- 9
notopdup$scale[notopdup$code==0873] <- 9
notopdup$scale[notopdup$code==0874] <- 10
notopdup$scale[notopdup$code==090] <- -2
notopdup$scale[notopdup$code==091] <- -2
notopdup$scale[notopdup$code==092] <- -2
notopdup$scale[notopdup$code==093] <- -2
notopdup$scale[notopdup$code==094] <- -2
notopdup$scale[notopdup$code==100] <- -5
notopdup$scale[notopdup$code==101] <- -5
notopdup$scale[notopdup$code==1013] <- -5#
notopdup$scale[notopdup$code==1014] <- -5#
notopdup$scale[notopdup$code==102] <- -5
notopdup$scale[notopdup$code==1031] <- -5#
notopdup$scale[notopdup$code==1033] <- -5#
notopdup$scale[notopdup$code==1041] <- -5
notopdup$scale[notopdup$code==1042] <- -5
notopdup$scale[notopdup$code==1043] <- -5
notopdup$scale[notopdup$code==1044] <- -5
notopdup$scale[notopdup$code==105] <- -5
notopdup$scale[notopdup$code==1053] <- -5#
notopdup$scale[notopdup$code==1056] <- -5#
notopdup$scale[notopdup$code==106] <- -5
notopdup$scale[notopdup$code==107] <- -5
notopdup$scale[notopdup$code==109] <- -5#
notopdup$scale[notopdup$code==110] <- -2
notopdup$scale[notopdup$code==111] <- -2
notopdup$scale[notopdup$code==112] <- -2
notopdup$scale[notopdup$code==1121] <- -2
notopdup$scale[notopdup$code==1122] <- -2
notopdup$scale[notopdup$code==1123] <- -2
notopdup$scale[notopdup$code==1124] <- -2
notopdup$scale[notopdup$code==1125] <- -2
notopdup$scale[notopdup$code==113] <- -2
notopdup$scale[notopdup$code==114] <- -2
notopdup$scale[notopdup$code==115] <- -2
notopdup$scale[notopdup$code==116] <- -2#
notopdup$scale[notopdup$code==120] <- -4
notopdup$scale[notopdup$code==121] <- -4
notopdup$scale[notopdup$code==1211] <- -4
notopdup$scale[notopdup$code==1212] <- -4
notopdup$scale[notopdup$code==1213] <- -4
notopdup$scale[notopdup$code==122] <- -4
notopdup$scale[notopdup$code==1222] <- -4#
notopdup$scale[notopdup$code==123] <- -4
notopdup$scale[notopdup$code==1231] <- -4
notopdup$scale[notopdup$code==1233] <- -4
notopdup$scale[notopdup$code==124] <- -5
notopdup$scale[notopdup$code==1241] <- -5#
notopdup$scale[notopdup$code==1243] <- -5#
notopdup$scale[notopdup$code==1244] <- -5#
notopdup$scale[notopdup$code==1246] <- -5#
notopdup$scale[notopdup$code==125] <- -5
notopdup$scale[notopdup$code==127] <- -5
notopdup$scale[notopdup$code==128] <- -5
notopdup$scale[notopdup$code==129] <- -5#
notopdup$scale[notopdup$code==130] <- -4.4
notopdup$scale[notopdup$code==131] <- -5.8
notopdup$scale[notopdup$code==1311] <- -5.8
notopdup$scale[notopdup$code==1312] <- -5.8
notopdup$scale[notopdup$code==1313] <- -5.8
notopdup$scale[notopdup$code==132] <- -5.8
notopdup$scale[notopdup$code==1322] <- -5.8
notopdup$scale[notopdup$code==133] <- -5.8
notopdup$scale[notopdup$code==134] <- -5.8
notopdup$scale[notopdup$code==136] <- -7
notopdup$scale[notopdup$code==137] <- -7
notopdup$scale[notopdup$code==138] <- -7
notopdup$scale[notopdup$code==1383] <- -7
notopdup$scale[notopdup$code==1384] <- -7
notopdup$scale[notopdup$code==139] <- -7
notopdup$scale[notopdup$code==140] <- -6.5
notopdup$scale[notopdup$code==141] <- -6.5
notopdup$scale[notopdup$code==1411] <- -6.5
notopdup$scale[notopdup$code==1412] <- -6.5
notopdup$scale[notopdup$code==1413] <- -6.5
notopdup$scale[notopdup$code==1414] <- -6.5
notopdup$scale[notopdup$code==142] <- -6.5
notopdup$scale[notopdup$code==143] <- -6.5
notopdup$scale[notopdup$code==144] <- -7.5
notopdup$scale[notopdup$code==145] <- -7.5
notopdup$scale[notopdup$code==150] <- -7.2
notopdup$scale[notopdup$code==151] <- -7.2
notopdup$scale[notopdup$code==152] <- -7.2
notopdup$scale[notopdup$code==153] <- -7.2
notopdup$scale[notopdup$code==154] <- -7.2
notopdup$scale[notopdup$code==160] <- -4
notopdup$scale[notopdup$code==161] <- -4
notopdup$scale[notopdup$code==162] <- -5.6
notopdup$scale[notopdup$code==1621] <- -5.6
notopdup$scale[notopdup$code==1622] <- -5.6
notopdup$scale[notopdup$code==163] <- -6.5
notopdup$scale[notopdup$code==164] <- -7
notopdup$scale[notopdup$code==166] <- -8
notopdup$scale[notopdup$code==1662] <- -8#
notopdup$scale[notopdup$code==170] <- -7
notopdup$scale[notopdup$code==171] <- -9.2
notopdup$scale[notopdup$code==1711] <- -9.2
notopdup$scale[notopdup$code==1712] <- -9.2
notopdup$scale[notopdup$code==172] <- -5
notopdup$scale[notopdup$code==1721] <- -5
notopdup$scale[notopdup$code==1722] <- -5
notopdup$scale[notopdup$code==1723] <- -5
notopdup$scale[notopdup$code==1724] <- -5
notopdup$scale[notopdup$code==173] <- -5
notopdup$scale[notopdup$code==174] <- -5
notopdup$scale[notopdup$code==175] <- -9
notopdup$scale[notopdup$code==180] <- -9
notopdup$scale[notopdup$code==181] <- -9
notopdup$scale[notopdup$code==182] <- -9.5
notopdup$scale[notopdup$code==1821] <- -9
notopdup$scale[notopdup$code==1822] <- -9
notopdup$scale[notopdup$code==1823] <- -10
notopdup$scale[notopdup$code==183] <- -10
notopdup$scale[notopdup$code==1831] <- -10
notopdup$scale[notopdup$code==1832] <- -10
notopdup$scale[notopdup$code==185] <- -8
notopdup$scale[notopdup$code==186] <- -10
notopdup$scale[notopdup$code==190] <- -10
notopdup$scale[notopdup$code==191] <- -9.5
notopdup$scale[notopdup$code==192] <- -9.5
notopdup$scale[notopdup$code==193] <- -10
notopdup$scale[notopdup$code==194] <- -10
notopdup$scale[notopdup$code==195] <- -10
notopdup$scale[notopdup$code==196] <- -9.5
notopdup$scale[notopdup$code==201] <- -9.5
notopdup$scale[notopdup$code==202] <- -10
notopdup$scale[notopdup$code==203] <- -10

sort(unique(notopdup$code)) # no NAs


################################################################

# assign verbal and material cooperation labels

notopdup$vcoop <- 0
notopdup$mcoop <- 0
notopdup$vcon <- 0
notopdup$mcon <- 0
notopdup$vcoop[notopdup$code >=1 & notopdup$code <=5] <- 1
notopdup$mcoop[notopdup$code >=6 & notopdup$code <=9] <- 1
notopdup$vcoop[notopdup$code >=10 & notopdup$code <=59] <- 1
notopdup$mcoop[notopdup$code >=60 & notopdup$code <=99] <- 1
notopdup$vcon[notopdup$code >=100 & notopdup$code <= 142] <- 1
notopdup$mcon[notopdup$code >=143 & notopdup$code <=204] <- 1
notopdup$vcoop[notopdup$code >=211 & notopdup$code <=356] <- 1
notopdup$mcoop[notopdup$code >=811 & notopdup$code <= 874] <- 1
notopdup$vcon[notopdup$code >= 1011 & notopdup$code <=1424] <- 1
notopdup$mcon[notopdup$code >=1431 & notopdup$code <=2042] <- 1

################################################################


head(top)
head(notop)
head(notopdup)


# create annual datasets

ann = data.frame(sort(unique(notop$year)))
colnames(ann) = "year"

ann[c("mconTop","mconNoTop","mconNoTopDup","scaleTop","scaleNoTop","scaleNoTopDup")] = NA
for(i in 1:nrow(ann)){
  ann$mconTop[i] = sum(na.omit(top$mcon[ann$year[i]==top$year]))
  ann$mconNoTop[i] = sum(na.omit(notop$mcon[ann$year[i]==notop$year]))
  ann$mconNoTopDup[i] = sum(na.omit(notopdup$mcon[ann$year[i]==notopdup$year]))
  ann$scaleTop[i] = mean(na.omit(top$scale[ann$year[i]==top$year]))
  ann$scaleNoTop[i] = mean(na.omit(notop$scale[ann$year[i]==top$year]))
  ann$scaleNoTopDup[i] = mean(na.omit(notopdup$scale[ann$year[i]==top$year]))
}

# 2 sided 10 year moving average
ma <- function(x,n=5){filter(x,rep(1/n,n), sides=2)}


plot(ann$year,ma(ann$mconNoTopDup),col="SlateGray",type="l",,xlab=NA,ylab="Annual no. material conflicts (10 year 2 sided moving average)")
points(ann$year,ma(ann$mconNoTop),col="coral4",type="l")
points(ann$year,ma(ann$mconTop),col="yellow",type="l")
legend("topleft",c("All articles","Less duplicates","Less sports, culture, obituaries"),col=c("SlateGray","coral4","yellow"),pch=16,bty="n")
dev.print(file=paste(figdir,"Sports1USA.pdf",sep=""), device=pdf); dev.off()


plot(ann$year,ma(ann$scaleNoTopDup),col="SlateGray",type="l",xlab=NA,ylab="Annual no. material conflicts (10 year 2 sided moving average)",ylim=c(-.5,3))
points(ann$year,ma(ann$scaleNoTop),col="coral4",type="l")
points(ann$year,ma(ann$scaleTop),col="yellow",type="l")
legend("topleft",c("All articles","Less duplicates","Less sports, culture, obituaries"),col=c("SlateGray","coral4","yellow"),pch=16,bty="n")
dev.print(file=paste(figdir,"Sports2USA.pdf",sep=""), device=pdf); dev.off()







